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Abstract 

A recent microstructure-based FEM model that couples crystal-based plasticity, the B2 <-> B19' phase transformation and anisotropic 
elasticity at the grain scale is calibrated to recent data for polycrystalline NiTi (49.9 at.% Ni). Inputs include anisotropic elastic prop- 
erties, texture and differential scanning calorimetry data, as well as a subset of recent isothermal deformation and load-biased thermal 
cycling data. The model is assessed against additional experimental data. Several experimental trends are captured - in particular, the 
transformation strain during thermal cycling monotonically increases and reaches a peak with increasing bias stress. This is achieved, 
in part, by modifying the martensite hardening matrix proposed by Patoor et al. [Patoor E, Eberhardt A, Berveiller M. J Phys IV 
1996;6:277]. Some experimental trends are underestimated - in particular, the ratcheting of macrostrain during thermal cycling. This 
may reflect a model limitation that transformation-plasticity coupling is captured on a coarse (grain) scale but not on a fine (martensitic 
plate) scale. 

Published by Elsevier Ltd. on behalf of Acta Materialia Inc. 
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1. Introduction 

An important challenge to shape memory alloy (SMA) 
models is to capture the actuation response [1-3]. A repre- 
sentative situation is shown in Fig. la, where the material is 
subjected to a constant tensile bias stress S b ias Upon heat- 
ing, the SMA contracts by a macroscopic transformation 
strain E T (feature F7). Other important features are the 
hysteretic width A(/ H (feature F6), widths A6 a V i (feature 
F8) and A0 m _ a for the austenite-to-martensite and mar- 
tensite-to-austenite transformations, respectively, and the 
open loop strain per thermal cycle, AE cyc[c These quantities 
typically depend on £ bias , composition, processing history 
and microstructural features such as single/polycrystalline 
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orientation/texture, precipitate morphology, dislocation 
substructure and internal stress (e.g. [4]). 

SMA models can be characterized by their microstruc- 
tural detail and phenomenological nature. For example, 
some track the volume fractions of stress-induced vs. ther- 
mally induced martensite via state variables [5,6] or adopt a 
functional dependence of transformation strain on the total 
martensite volume fraction and stress [7]. Others include a 
phenomenological description of transformation-induced 
plasticity [8-11]. Variant-level models (e.g. [12-20]) track 
the volume fraction of individual martensite variants, 
thereby capturing effects of texture and crystal orientation. 
One version includes a recursive algorithm method that 
adjusts interactions between competing tetradomains in 
order to reproduce neutron diffraction data [21], The 
method captures the evolution of martensite volume frac- 
tion and local stress for various (hkl) sets of domains [22]. 
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Fig. 1. Axial macrostrain E vs. temperature 0 during thermal cycling with an axial bias stress E bias : (a) experimental data for the 49.9 at.% Ni-Ti alloy and 
(b) calibrated model results. The inserts show the transformation strain E T vs. E b ; as , where E T is defined by feature F7 in (a). The calibrated model 
parameters are presented in Table 1. 


A particular challenge to simulate actuator performance 
is to capture the initial, gradual increase in E T with £ bias , 
the attainment of a peak E T at moderate E b j as and the 
decrease in E T at larger £ bias (Fig. la, inset). The gradual 
increase in E T tests model assumptions about martensite 
interactions. These can be formulated between individual 
twinned habit plane variants (hpv) [14-17,23], hpv groups 
[24] and lattice correspondence variants (lev) [25], or 
groups thereof [21,22]. An implementation of the first type 
by Patoor et al. [23] has been used (e.g. [15-17]), but the 
predicted E T increases too abruptly with X bias [25]. 

The peak in E T with S bias tests model assumptions con- 
cerning plasticity. Most do not incorporate it at all [12-16], 
including the aforementioned recursive model [21,22]. 
Some recent efforts have included martensite plasticity, 
either in J2-based [18,26] or crystal-based [17] forms. How- 
ever, austenite rather than martensite plasticity has been 
observed in recent transmission electron microscope stud- 
ies of pseudoelastically deformed, solutionized 50.7 at.% 
Ni-Ti single crystal pillars [27] and cold-worked 50.8 
at.% Ni-Ti wires [28], as well as thermally cycled, single- 
crystal 50.4 at.% Ni-Ti [29]. In particular, a detailed anal- 
ysis reveals that transformation-induced dislocations have 
Burgers vectors and slip planes with rationale directions 
in the austenite basis rather than the martensite basis 
[27]. These observations follow earlier work documenting 
the formation of planar dislocation arrays during load- 
biased thermal cycling of 50.1 and 50.4 at.% Ni-Ti [30]. 
Plastic deformation has also been proposed to explain pro- 
gressive widening of the { 1 1 0} B2 peak during load-biased 
thermal cycling [31]. 

Consequently, a principal contribution of this paper is 
to determine whether a crystal-based modeling approach 


that incorporates austenite plasticity can capture the key 
thermal cycling features shown in Fig. la as well as other 
data. This recent finite element (FE) approach [32] couples 
a crystal-based description of the evolution of martensite 
hpv by Thamburaja and Anand [13] with a crystal-based 
description of austenite plasticity based on Peirce et al. 
[33]. Thus, it is able to assess the Patoor et al. [23] formal- 
ism for hpv-hpv interaction as well as hpv-crystal plastic- 
ity interaction. Other unique features are inclusion of 
rigorous formalisms for anisotropic elastic deformation in 
the austenite and martensite phases, as well as predictions 
of the orientation of individual Bain strain variants of mar- 
tensite. Recently reported elastic constants for monoclinic 
martensite [34] are included, as well as the appropriate 
invariant-plane condition between austenite and martensite 
- a feature not satisfied in a recent self-consistent approach 
to incorporate anisotropic elasticity [35]. 

The structure of the manuscript is as follows: Section 2 
describes experimental data for the polycrystalline 49.9 
at.% Ni-Ti material system, performed in support of the 
modeling activity. Section 3 describes the FE model. Sec- 
tion 4 presents the process for calibrating the model param- 
eters to a subset of experimental data, and Section 5 
assesses predictions of the calibrated model to additional 
data. Section 6 provides an overall assessment of model 
capabilities and conclusions. 

2. Experimental characterization 

2.1. Material system and history 

The binary 49.9 at.% (55 wt.%) Ni-Ti alloy used in this 
study is described elsewhere [2], It was produced by Special 
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Metals, New Hartford, NY and supplied as 10 mm diame- 
ter rods in a hot-rolled/hot-drawn and hot-straightened 
condition. The specific thermo-mechanical processing his- 
tory is proprietary. The alloy is single phase, with a dynam- 
ically recrystallized and equiaxed grain structure with 
~40 pm average grain size. 

All experimental measurements on as-machined speci- 
mens were preceded by two no-load thermal cycles from 
30 to 200 °C at 20 °C min 1 . The nominal purpose is to 
reduce internal stress generated by material processing 
and machining. 

2.2. Texture analysis 

The as-received alloy had a very weak recrystallization 
texture. Fig. 2a shows the pole figures in the austenite 
phase at 200 °C. They were obtained using the high-pres- 
sure preferred orientation (HIPPO) neutron diffractometer 
at Los Alamos National Laboratory (LANL) [36]. 

2.3. Differential scanning calorimetry 

Differential scanning calorimetry (DSC) tests were per- 
formed to measure the heat of transformation, £> DSC , the 
martensite start and finish temperatures, 0 M s and 0 MF , 
and the austenite start and finish temperatures, 0 AS and 
Oaf [2]. Key results were determined from a« = 7 DSC 
sampling size and are summarized in Section 4.2. 



(a) Experiments (b) simulations 

Fig. 2. Pole figures for as-received austenitic polycrystalline 49.9 at.% Ni- 
Ti at 200 °C from (a) neutron diffraction data using HIPPO and (b) 
simulations based on SMARTS diffraction data [2] for the same material 
at 200 °C, using POPLA software [46] and assuming axisymmetry. A 
stronger (111) and weaker (100) texture are apparent. 


2.4. Pseudoelastic and plastic strains above 0 AF 

Fig. 3a shows the uniaxial tension stress-strain 
response at 0 o =13O°C for maximum imposed axial 
strains £ mas = 0.02 and 0.04, and also at 215 °C for 
fs max = 0.04. Stress-induced martensite forms during 
deformation above 0 AF and some is retained upon 
unloading. A small axial strain rate of 10 4 s 1 was used 
to help ensure isothermal conditions. After F max was 
attained, the sample was unloaded isothermally and then 
heated to 600 °C under no load. During heating, the 
macrostrain E decreased due to recovery (transformation) 
of stabilized martensite to austenite. It was then cooled 
back to 0 O Recovery was confirmed by noting that the 
macrostrain-temperature (E-8) response during initial 
cooling is consistent with the thermal expansion coeffi- 
cient of austenite [37]. The strain £p OSt -heat (feature FI 3, 
Fig. 3a) at 0 O is expected to be due primarily to plastic 
deformation. 

This NiTi alloy does not exhibit a strong, “flag-like” 
pseudoelastic response, as typically observed above 0 A f 
in Ni-rich NiTi alloys [12-15,28,29,38M2] or cold-worked 
and annealed Ti-rich NiTi alloys [43]. The more idealized 
pseudoelastic response, with clear loading and unloading 
plateaus and large recoverable strains, is obtained in micro- 
structures that suppress plastic deformation [44] yet permit 
stress-induced transformation. In practice, this is achieved 
through grain refinement [28,42], cold work with partial 
recovery [43] or precipitation strengthening of Ni-rich NiTi 
alloys through suitable heat treatment [38-41]. Such strat- 
egies were not pursued here since the goal was to study the 
simultaneous effects of plasticity and transformation. 
Indeed, the features in Fig. 3a (and Fig. la) show evidence 
of both plasticity and pseudoelastic transformation, mak- 
ing it a suitable system to test the polycrystalline 
simulations. 

2.5. Load-biased thermal cycling 

Fig. la shows the E-8 response of the NiTi alloy at dif- 
ferent values of uniaxial bias stress (X bias ). The stress level 
was achieved by imposing a 10 4 s ' axial strain rate at 
0 m in = 30 °C, where martensite is stable, until the desired 
stress was reached at constant temperature. The controller 
was switched into load control and the stress was held con- 
stant. This was followed by thermal cycling up to 
0 ma x = 200 °C and back to 0 min = 30 °C at constant X bias 
The second thermal cycle responses are shown in Fig. la, 
as well as the transformation strain E T vs. Z bias (inset). 

2.6. In situ neutron diffraction 

Neutron diffraction spectra were obtained during stress- 
biased thermal cycling, using the protocol in Section 2.5. 
These measurements were obtained in “time-of-flight” 
mode using the Spectrometer for MAterials Research at 
Temperature and Stress (SMARTS) facility at LANL. 
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Fig. 3. Axial macrostress X vs. macrostrain E at Oo = 130 vs. 215 °C showing (a) experimental data for the 49.9 at.% Ni-Ti alloy and (b) calibrated model 
results, ^unload and 2? post . heat are the respective macrostrains after unloading and after a subsequent 600 °C heat treatment of the unloaded sample. The 
calibrated model parameters are presented in Table 1. 


Two detector banks furnished the diffracted intensity vs. d- 
spacing from two groups of crystallographic planes: those 
parallel and those perpendicular to the loading axis. The 
data give the evolution of martensite and austenite (volume 
fraction and orientation) during thermal cycling. Addi- 
tional details are provided in Ref. [2] and the references 
therein. 

3. Polycrystalline simulations and constitutive relation 

3.1. Polycrystalline simulations 

A brief description of the finite element representation 
of the polycrystal and boundary conditions is provided, 
with emphasis on texture specification and new model 
extensions as noted in Section 1 . Additional details are pro- 
vided in Ref. [32]. 

3.1.1. Discretization at the grain scale and texture 
specification 

The polycrystal was modeled by a 7 x 7 x 7 cube assem- 
bly of 8-node three-dimensional brick elements (C3D8 for 
isothermal analyses and C3D8T for thermo-mechanical 
analyses) using the commercial FE software ABAQUS 
[45]. Each element represented a grain with an assigned 
crystal orientation of the austenite (B2) phase. These orien- 
tations were based on a simulated texture for the as- 
received NiTi (Fig. 2b), obtained by processing SMARTS 
data [2] with the texture analysis software POPLA [46], 
assuming an axisymmetric texture. This assumption and 
the stronger (111) and relatively weaker (10 0) textures 
present in the SMARTS data compare favorably with sub- 
sequent ELIPPO data shown in Fig. 2a (Section 2.2). The 
number of grains, /V grain = 343, in the simulations is suffi- 
cient since the macroscopic response does not change sig- 
nificantly for Agrain > 343. 


3.1.2. Predeformation E p ( pre ) in the austenitic state 

Some internal stress is expected in the samples prior to 
testing, even though the material was hot-worked and 
two no-load or stress-free thermal cycles were imposed 
before testing. The simulations introduced internal stress 
via predeformation in the austenitic state. This was 
achieved by heating to 300 °C, at which temperature aus- 
tenite is stable, imposing a macroscopic plastic prestrain 
Fp(pre) at 10 4 s ', unloading, then cooling to the desired 
test temperature 0 o E p(pre) was varied to match the experi- 
mental data (Section 4.3). 

3.1.3. Isothermal deformation testing and post-lieating 

The polycrystal was strained to E m;lx at some constant 

temperature 0 o > 0 AF and then unloaded. This was 
achieved via prescribed normal displacement rates and zero 
shear tractions on the top and bottom surfaces, equivalent 
to an axial strain rate of ±10 _4 s~ 1 . The transverse faces 
were traction free. The macrostrain Funioad after unloading 
was partitioned into plastic and transformation contribu- 
tions by integration of the local plastic and transformation 
strains over the polycrystal. The simulations assumed spa- 
tially uniform, isothermal conditions. This is considered 
reasonable since the small imposed strain rate minimizes 
local heating and heat transfer rates associated with the 
austenite-martensite phase transformation. Fig. 3b shows 
results using the calibrated model. 

3.1.4. Stress-biased thermal cycling 

The polycrystal was heated to the maximum cycling 
temperature 0 max =190 °C, then the bias stress E bias was 
imposed using an axial strain rate of 10 4 s '. The sample 
was then cycled between 0 min = 30 °C and 0 max at 
±0.08 °C s '. The simulations assumed a spatially uniform 
temperature, which is reasonable, given the small heating/ 
cooling rate. The simulations began the load application 
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and cycling at 0 max , compared to 0 min for the experiments. 
However, previous experiments [3] and the current results 
show this does not affect key parameters such as transfor- 
mation strain. Fig. lb shows results using the calibrated 
model. 


3.1.5. In situ neutron diffraction simulations 

Three post-processing steps were performed to compare 
simulations and experiments. First, the predicted volume 
fractions of martensite habit plane variants were converted 
into volume fractions of martensite correspondence vari- 
ants and their orientation using the Crystallographic The- 
ory of Martensite (CTM) [47] - in particular, Eq. (17) of 
Ref. [48] was solved. Second, the martensite volume frac- 
tions V]v|| and v w± were determined, where N denotes the 
plane normal of interest (e.g. [1 0 0] or [0 1 1]) and || and 
1 denote whether N is parallel or perpendicular to the 
loading axis. Third, the diffracted intensity 7(N) of plane 
N was normalized by the intensity / re f{N) from a reference 
test case so that 


/( N) ^ v N 
7 re /(N) v Nref 


( 1 ) 


Load-biased thermal cycling under a bias stress of 
100 MPa was chosen as the reference. 


3.2. Constitutive relation for austenite-martensite 
aggregates 

The finite element formulation requires a relation 
between an imposed deformation gradient on an austen- 
ite-martensite aggregate and the corresponding increments 
in average stress, martensite volume fraction and slip sys- 
tem activity. For simulation of entire polycrystals, it is 
impractical to retain detailed spatial information concern- 
ing individual martensite variants; instead, the aggregate 
is described by the average volume fraction of martensite 
habit plane variants (plates). The averaging process can 
be viewed as an isostress approach in that elastic, plastic 
and transformation processes are driven by a common 
aggregate Cauchy stress T and they collectively contribute 
to the average deformation gradient of the aggregate. 
These features and improvements to the original formula- 
tion [32] are highlighted below. 


martensite phase, following Ref. [50]. First, the symmetric 
Piola-Kirchoff stress T* and Cauchy stress T are related by 


T = J L.. F e T*F eX where T* 
det(F ) 

= C[E e — A th (0 — 0 O )] 


di// 

ar 


( 2 ) 


T r denotes the Helmholtz free energy of the austenite-mar- 
tensite aggregate as defined in Eq. (8) of Ref. [32], and the 
average elastic strain is E e = !/ 2 (F eT F e - I). Eq. (2) requires 
the aggregate thermal expansion coefficient A th : 

Ath = VAAth-A + Ath— M (3) 

and effective elastic moduli C: 


n t 

C = v A C A + ^2 v,Cmw where 

7=1 ( 4 ) 

Cm(i) = V/,/Cm(<,;) + v ( jCm(<j) 

Austenite has cubic symmetry so that Ca has three inde- 
pendent components. Each habit plate variant t has a local 
volume fraction v, and elastic modulus C M(n , given by the 
average moduli C M(M) , C M (^ and volume fractions v tJ , v tJ 
(=1 — v t j) of the monoclinic (B19') correspondence vari- 
ants i and j within plate t. For practical purposes, the com- 
ponents are referred to the austenite crystal basis. This 
requires the rotation operators Q t i and Q, ; to relate a 
direction in the martensite crystal basis of variants i and 
j, respectively, to the crystal basis of austenite from which 
plate t formed. They are obtained by solving Eq. (17) from 
the CTM [48]: 

Q I,i = R/R/Rlc; and Q tj = R/R|C, (5) 

where R L c-/ and R|,c-/ are the lattice correspondence rota- 
tions between martensite variants i and j and austenite, 
respectively, R (/ is the rotation required to form a twin 
plane between variants i and j, and R (/ is the rotation re- 
quired to form a habit plane with the austenite [48,50]. 
Although the thermal expansion coefficient for martensite 
plates can be constructed in a similar fashion, the isotropic 
form A th _ M = a M I is adopted with the supposition that the 
primary effect of anisotropy stems from the elastic moduli. 
This assumption is discussed further in Section 5.2.4 and 
justified by the small magnitude of thermal strain 
(~0.1%) during thermal cycling [37]. 


3.2.1. Partitioning elastic and inelastic deformation 
Following Lee [49], the deformation gradient F = F e 

F mel , so that it is multiplicatively decomposed into a 
thermo-elastic part F e and an inelastic part F mel . Sec- 
tion 3.2.2 describes how F e is related to the Cauchy stress 
T and Sections 3. 2. 3-5 describe how the time derivative 
of inelastic deformation F mel is related to T. 

3.2.2. Elastic constitutive relation 

This section highlights modification of the existing 
theory to incorporate the full elastic anisotropy of the 


3.2.3. Average inelastic velocity gradient for the aggregate 
The time derivative of inelastic deformation is 
pme! _ pmei pnei^ w j iere q lc inelastic velocity gradient of 
the austenite-martensite aggregate is approximated by a 
rule of mixtures [32]: 

N§ Nj 

L mel = v A 2 A( b f P ® m f P ) + 2 v, (b^ rans ® m ] 1 ™) ; 

S=1 ' =1 ( 6 ) 

N t 

V A = 1 - v m and v M = ^ v, 

7=1 
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The first term is the contribution from plastic deformation 
of a volume fraction v A of austenite having N s slip systems 
with Burgers vectors bf p , slip plane normals mf p and rates 
y s of slip activities. The second term is the contribution 
from the rate of increase v, in volume fraction of martensite 
plate type t (=1 -N T ) with habit plane normal m| rans and 
transformation displacement bj rans . The second equation 
acknowledges that the austenite plus martensite volume 
fractions equals 1, and the third equation indicates that 
the total martensite volume fraction is the sum of the plate 
volume fractions. y s and v, are related to the Piola-Kirchoff 
stress T* in Sections 3.2.4 and 3.2.5, respectively. 


3.2.4. Plastic constitutive relation 

The slip rate % in Eq. (6) is specified according to Peirce 
et al. [33]: 


ys = y 0 


(1/m) 


sign(T ! );T J = bf p - (F eT F e r) -mf p 


( 7 ) 


where t s . is the resolved shear stress on slip system s, y 0 is a 
reference shear rate, m is the rate sensitivity and g s is the 
slip system hardness that evolves to a saturation value g sat 
according to 


Ns 

gs = XA(2+ ( i - Q) s *r)\y r \ 

r= 1 


and h r = /z 0 



( 8 ) 


The self and latent hardening rates are h r and Qh r , 
respectively. 


3.2.5. Transformation constitutive relation 

The rate of change v, is determined by ensuring that the 
unit volume driving force f, to form plate type t is bounded 
as —f c ^ f < / c , where f c is a critical value and 

F f (v„ T, 9) = bj aas • (F eT F e r ) • - g (9 - 0 T ) 

Nj 

~Y^Kv u (9) 

u= 1 

This formalism, proposed by Thamburaja and Anand [13], 
specifies no change in v, when f, satisfies the inequalities 
preceding Eq. (9). Should f > / c , v, is increased so as to ren- 
der f t =f c and if f < —f c , v t is decreased so as to render 
ft — fc- Thus, it is similar to a plastic yield surface.The driv- 
ing force ft has three terms. The first is the mechanical work 
of the aggregate stress acting through the transformation 
strain induced by the formation of plate type t. The second 
is the chemical energy associated with the latent heat of 
transformation per unit reference volume, 1 T . The third is 
the interaction between plates in the aggregate. It repre- 
sents, for example, the additional contribution to f due 
to the stress field of nearby martensite plates u; such de- 
tailed interactions are not captured by the aggregate stress 
in the first term. 

4. Calibration of material parameters 

Table 1 shows the material parameters used to simulate 
the 49.9 at.% (55 wt.%) Ni-Ti SMA. They are adopted 


Table 1 

Calibrated properties for hot rolled/hot drawn polycrystalline 49.9 at.% Ni-Ti (55wt% Ni). 


Austenite structure 

(see Section 3.1.1) 


Property 

Value 

Calibration feature 

Texture 

Pole figure (see Fig. 2b) 

FI: Fig 2a 

Elastic and thermal, Ai-thermai = (Ca, C m . ^th-A, ^th-ivil (see Section 4.1) 


C A (GPa) 

C„ = 183 C n = 146 C44 = 46 

F2: DFT calculations [34] 

C M (GPa) 

C,, = 249 Cn = 129 C13 = 107 C15 = 15 C 22 = 245 C 23 = 125 
C 25 = -3 C33 = 212 C35 = -1 C44 = 87 C 46 = -4 C 55 = 66 
C 66 = 86 Hill Averages: E = 180, G = 69 

F2: DFT calculations [34] 

«a(K ') 

1 1 x 10~ 6 

F3: Neutron diffraction [37] 

«m(K ') 

6.6 x 10~ 6 

F4: Neutron diffraction [37] 

Phase transformation, P lrins = (0 T , 2 X , f c and h, u } (see Section 4.2) 


MK) 

353.7 

F5: DSC data 

2 t (MJ m 3 ) 

140 initial guess: 165.6 

FI 1-13: Fig. 3a F5: DSC data 

fc (MPa) 

0.05U t (=7.14) 

F6: Fig. la and Eq. (10): 

h, u (MPa) 

hcom — 0; /line — Ca( 44 )/ 12,000; /?self — Ca(44)/400 

F7-8: Fig. la. F9: Fig. 3a and (11) See also Table 2 [15] 

Austenite plasticity, P p i as tio = lyo, g sa t, g?, K- Q - a, -Epipre)} (see Section 4.3) 


7o (s ) 

0.002 

F10: literature [32] 

m 

0.02 

F10: literature [32] 

Q 

1.4 

F10: literature [32] 

a 

0.125 

F10: literature [32] 

g sa , (MPa) 

900 

F10: literature [32] 

g“ (MPa) 

272 

FI 1-14: Fig. 3a 

h 0 (MPa) 

50 

FI 1-14: Fig. 3a 

^p(pre) (%) 

-0.7 

FI 1-14: Fig. 3a 
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Fig. 4. (a) Axial macrostrain E vs. temperature 0 during thermal cycling with an axial bias stress Xbias = 50 MPa and (b) axial macrostress £ vs. 
macrostrain E at 0 o = 130 °C for the polycrystalline 49.9 at.% Ni-Ti alloy. The calibrated model parameters are in Table 1. Other results in (a) use the 
calibrated parameters except with the interaction matrix h,„ of Patoor et al. [23] or with E p(pre) = 0. The £ p(pre| = 0 cases in (b) are nearly coincident and 
use (g°, h 0 ) = (235 MPa, 500 MPa) vs. (250 MPa, 50 MPa) (the former is slightly higher). 


from published literature values or calibrated to DSC data 
and select isothermal deformation and thermal cycling tests 
performed specifically to support this modeling effort. 
These various “features” and the calibration process are 
described in Sections 4. 1-4.3. 

4.1. Elastic-thermal, P e i. thermu i- {C A , C M , A th _ A , A th . M ) 

The anisotropic elastic stiffness C A for the cubic (B2) 
austenite phase and C M for the monoclinic ( B 19') martens- 
ite phase are adopted from the density functional theory 
calculations (feature F2, Table 1) of Hatcher et al. [34], Iso- 
tropic thermal expansion properties are assumed with a A 
and a M taken from the neutron diffraction experiments 
(feature F3, Table 1) of Qiu et al. [37]. Use of recently com- 
puted elastic moduli corrects frequent assumptions in the 
literature that C M ~ V£C A (e.g. [12,13]) or C M < C A (e.g. 
[7,9,15-18], These effects, including differences in moduli 
from Hatcher et al. [34] vs. Wagner et al. [51], are discussed 
in Section 5.2.4. 

4.2. Transformation. P, rcms = {0 T , ). T , f c and 

This study models /V T = 24 type-II twinned martensite 
plate types, like many micromechanics-based models 
[13,15-17,26], although there are 192 theoretically possible 
habit plane variants (hpv) for B2 — ► B19' transformation 
[47]. The transformation parameters Arans are obtained 
from a combination of DSC (no load), isothermal deforma- 
tion tests and thermal cycling data at Z hias = 50 MPa. The 
DSC data (averaged over seven samples) gives 0 M f> d MS , 
0 AS and 0 AF = 46, 71, 86 and 109 ± 2 °C, respectively, for 
this 49 at.% Ni-Ti alloy. These values (feature F4) furnish 
the transformation temperature 0 T = (0 M s + 0 AS )/ 
2 = 353.7 K (standard deviation = 2 K). The calibrated 
simulations use a latent heat of transformation 


A t = 140 MJ nT 3 , as detailed in Section 4.3. Initially, an 
average of Qusc (feature F5) for the forward A — ► M 
and reverse M — > A heats of transformation is used, 
yielding /. T « Qdsc = 165.6 MJ nT 3 (standard deviation = 
4.7 MJ m 3 ). This is viewed as an initial guess since £>dsc 
includes contributions from elastic or defect energies associ- 
ated with the transformation. 

The critical driving force for transformation is obtained 
from 


f c _ 0 M ^ A - 0 A ^ M _ A0 h 

Aj 2 6j 2 0j 


( 10 ) 


Eq. (10) is obtained by writing the forward (/ A ^ M = f c ) 
and backward (/] M_ * A — —fc) critical conditions in extended 
form using Eq. (9) and then taking the difference between 
them. Eq. (10) furnishes f c /A T = 0.051, based on an average 
hysteretic width (feature F6) A0 H = 36.3 K (standard devi- 
ation = 1.15 K) for thermal cycling at £ bias = 50 MPa (see 
Fig. la). This small £ bias minimizes plasticity; a similar fj 
i.y value occurs with E bias = 0 data. Formally, A 0 H is an 
average width over the range £ = 0. 1-0.9, where C is the 
fraction of transformation strain. 

The structure of the martensite plate interaction matrix 
h tu in Eq. (9) is obtained by calibrating the simulations to 
give E t « 1% (feature F7) and A0 a ^ m « 25 °C (feature 
F8) at E b i as = 50 MPa. A0 a ~ m is the decrease [Qmf—Oms) 
in temperature to complete the A — > M transformation. 
The outcome is that h, u must have three independent 
values: 

r Kck(= Ca ( 44)/400) if t = u 

h, u = < h com (= 0) if t A u and det(e T (o - e T („)) = 0 

[ h a (= Ca(44)/12000) if t A u and det(er(<) — £t ( q) ¥= 0 

( 11 ) 


This structure is a modihcation to the popular form 
proposed by Patoor et al. [23], where /j se i f = h com = 
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C A (44)/3000 MPa and h inc = C A(44) /7 50 MPa are used for 
NiTi [15-17]. The rationale for the modification is that sim- 
ulations employing the Patoor et al. [23] structure (curve 

Patoor), Fig- 4a) give E T ~ 4% for E bias = 50 MPa, and 
thus over-predict the experimental values, regardless of 
the assumed C A ( 44 , value. Martensite plates with maximal 
transformation strain along the loading axis dominate over 
other plate types because the stress work (term 1 in /), Eq. 
(9)) initially renders the driving force /+ largest for this fa- 
vored (+) plate type. As long as h com < h mc , the favored (+) 
plates experience less hardening than non-favored (-) 
plates, because -df+/dv + < df-/dv + . The modification 
sets /; se if » /; CO m, h i nc to suppress formation of a single plate 
type. 

The calibration is guided by application of Eq. (9) to an 
idealized case of two populations of plate types: a favored 
(+) type that renders a transformation strain E T(max) if 
v+ — > 1 and a non-favored (-) collection of self-accommo- 
dating plate types that renders E T ~ 0 if v_ — > 1. Taking 
the difference between the critical force conditions /+ = f c 
and /_ = f c furnishes 


sv s = 

bias 


A h+- 
2 




( 12 ) 


where S t Z b ias is the aggregate stress work to form favored 

plates and A A + _ ■ h ++ — /z_ + and A /j_ + = h — /?+_ are 

differences in interaction elements. Thus, increases in the 
differences (A se if-A com ) and (A S eir-A inc ) increase the required 
S i I h i;is for favored plates, decreasing E v . Eq. (9) also 
furnishes 


. 0 MS — Amf 



max MS S+ 


A h 


++ 


A/z h + A/z_ + 


(5+ Z +A A_+) + A+_ 

bias 


and 

(13) 


where h++ = h++ — h + _. The second equation of Eqs. (13) 
suggests that A self can be determined from the difference, 
S ma x - Z MS (feature F9), from isothermal deformation 
tests (Fig. 3a). Further, Eq. (12) and the first equation of 
Eqs. (13) suggest that differences such as h sel[ -h com and 
A S eif-A inc can be determined by capturing Features F7 
and F8 for thermal cycling at low bias stress (Fig. la). 
The final results reported in Eq. (11) and Table 1 are ob- 
tained as a best fit. 


4.3. Austenite plasticity, P p i astic = {y 0 ,m, g° s , g sat , h 0 , Q, a, 
Ep(pre ) } an d X T update 

These parameters are obtained from literature values 
and additional isothermal deformation data. In particular, 
6{110}/(10 0) and 6{0 0 1}/(10 0) austenitic slip sys- 
tems are considered [52]. The rate sensitivity m = 0.02, ref- 
erence strain rate y 0 = 10“ 3 s ' (E ^ UP 4 s '), g=1.4, 
a = 0.125 and g sat = 900 MPa (see Eqs. (7) and (8)) are 
adopted from prior calibration of the model to single crys- 
tal pseudoelastic compression tests (feature F10), using 


loading orientations that enhance plasticity [32]. In reality, 
there is limited data to determine these parameters accu- 
rately, but the simulations are less sensitive to them. The 
remaining plasticity parameters are the threshold resolved 
shear stress g° s for plastic flow, initial strain hardening h 0 
and plastic pre-deformation 7s p(pre) . The resulting calibra- 
tion gives a best fit of g° s = 272 MPa, h 0 = 50 MPa, E p(pre) 
= —0.7% and 1 T = 140 MJ m 3 (updated from an initial 
guess 165.6 MJ m ’). 

More specifically, the calibration procedure matches 
E \ i s ? E max and Epost-heat 0.0085 for E max 0.02 (Fea- 
tures Fll-13, Fig. 4b). The updated, smaller a t decreases 
Xyis and S max to better match the data, but g° s and h 0 are 
indeterminant. For example, the same curve ( E p{ pre) = 0, 
Fig. 4b) is obtained for (g®, h 0 ) = (250 MPa, 50 MPa) vs. 
(235 MPa, 500 MPa). The first case is selected since it bet- 
ter matches E post _ heat = 0.0204 (feature F14, Fig. 3a). The 
second case under predicts the value. 

A final issue is that curve /: p(prc) = 0 (Fig. 4b) over pre- 
dicts S M s (feature Fll). Prior work [32] demonstrates that 
compressive prestraining decreases S MS and increases 
hardening. The “calibrated model” result (Fig. 4b) adopts 
E p( pre) = -0.7% and (g° s , h 0 ) = ( 272 MPa, 50 MPa) as a 
best match. A larger compressive prestrain does not 
improve the match. This pre-deformation produces a 
multi-axial and spatially varying residual stress distribution 
with stress in the loading direction varying from —345 MPa 
to 272 MPa. These magnitudes are a substantial fraction of 
Z MS (feature Fll) in Fig. 3a. 

5. Model assessment 

5.7. Isothermal deformation response 

The calibrated model captures the monotonic loading 
paths for different E max and 6. In particular, £ max for 
6 = 130 °C (Pt A, Fig. 3b) agrees with experiments. Also, 
both experiments and simulations show an upward shift 
A£ ~ 100 MPa when 0 increases to 215 °C (Fig. 3b). If 
plasticity is “turned off” in the simulations, an excessive 
shift AS = 627 MPa occurs. This is consistent with a 
Clausius-Clapeyron slope dS/d 0 = 6.2 MPa K 1 for 
favorably oriented martensite plates [15,53]. The simula- 
tions show that plasticity intervenes to lower the shift; it 
is driven by complex, multi-axial stress states inherent in 
polycrystals. 

The discrepancy between the experimental and simu- 
lated unloading paths suggests that the reverse (M — ► A) 
transformation is complex. The experimental unloading 
paths are relatively linear (Fig. 3a), with most of the 
stress-induced martensite recovered during subsequent 
heating rather than during unloading. The simulated 
unloading path (6 = 130 °C, Fig. 3b) is the opposite, with 
most of the stress-induced martensite recovered during 
unloading rather than heating. A possible interpretation 
is that the simulations underestimate f c for the reverse 
transformation. This might be due to coupling between 
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plasticity and transformation at fine (variant) scales - a 
phenomenon not captured explicitly by the aggregate 
formulation. 

5.2. Thermal cycling response 

Overall, several experimental trends are captured, 
including the textured nature of martensite, but again there 
are quantitative discrepancies. 

5.2.1. Effect of bias stress and prestrain on transformation 
strain and critical temperatures 

Consistent with experiments, the predicted transforma- 
tion strain E T increases with E bias , as do the critical temper- 
atures 9 as , 9 af , 9 ms and 0 M f> and the temperature 
difference 0 M s — 6 mf (Fig. 1). Fig. 5 shows that even the 
two-way response for thermal cycling at I bias = 0 is cap- 
tured. This demonstrates that a plastic predeformation 
Fp(pre) — -0.7% creates a positive E T , by inducing an inter- 
nal stress state that biases martensite plate formation. Con- 
versely, a negative E T is generated for /: p(prc| > 0. It is 
interesting that a nonzero E T occurs, even though the aver- 
age internal stress is zero. 

Fig. 6 shows that the simulations capture a peak in ll\ 
with E bias . This is due primarily to the new proposed struc- 
ture for h tu (Eq. (11)) and the coupling of crystal plasticity 
and transformation (Eq. (6)). These features cause E T to 
increase gradually with £ bias due to formation of the 
favored (+) martensite. For E bias > 350 MPa, the model 
predicts a decrease in E T with X bias due to the onset of open 
loop strain, AE cycle (see Fig. 8). The predicted effect of 
Fp(pre) = -0.7% vs. 0 on the E i— E bias response in Fig. 6 is 
negligible. For comparison, simulations using the Patoor 
et al. [23] structure for h tu substantially overestimate E T 
at small E bias (curve h tu( p at0 or), Fig. 6) and predict E T to 
be relatively independent of Z bias , contrary to these and 
previous experimental results [1-3,30]. 



Temperature 0 (°C) 

Fig. 5. Axial macrostrain E vs. temperature 0 during no-load thermal 
cycling of the polycrystalline 49.9 at.% Ni-Ti alloy, showing a two-way 
effect. The calibrated model parameters are presented in Table 1. 



Fig. 6. Axial transformation strain E T vs. bias stress Ebias during thermal 
cycling between 0 m ; n = 30 °C and 0 max = 165 °C (lower experimental 
curve) vs. 200 °C (upper experimental curve) for the polycrystalline 49.9 
at.% Ni-Ti alloy. The calibrated model parameters are presented in 
Table 1. Also shown is the calibrated model result using the interaction 
matrix h lu of Patoor et al. [23], The model results are insensitive to 
0 max — 165 VS. 200 °C. 


5.2.2. Plastic strain enhancement due to phase 
transformation 

Simulations of thermal cycling at large £ bias (400 MPa) 
reveal that transformation enhances the macroscopic plas- 
tic strain E p . Fig. 7 shows rather modest increases in E p 
with thermal cycling for the “plasticity only” case - where 
transformation is turned off. In contrast, E p is larger in the 
“calibrated model” case, where both transformation and 
plasticity are present. The largest increases in E p occur dur- 
ing the early stages of cooling and latter stages of heating, 
when v M is small (see regions 1 and 2, Fig. 7). The phenom- 
enon is reminiscent of ratcheting of E p during stress-biased 
thermal cycling of composites with a large thermal strain 
mismatch between constituents [54]. These predictions are 
consistent with prior work [32] showing that plasticity aids 



Fig. 7. Predicted axial plastic macrostrain E p vs. temperature 0 during 
thermal cycling with an axial bias stress Z bias = 400 MPa. The calibrated 
model parameters are presented in Table 1. The "plasticity only” case uses 
the calibrated model without any phase transformations. 
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Temperature 0 (°C) 

Fig. 8. Axial macrostrain E vs. temperature 0 during thermal cycling with 
an axial bias stress E bias = 400 MPa, showing predictions of the calibrated 
model, which uses elastic constants C M and Ca from Ref. [34], Also shown 
are the calibrated model predictions with C M ~ Zi C A , C M and C A from 
Ref. [51], and C M = C M(Isotropic) . The last uses an isotropic matrix with 
Young’s modulus E and elastic shear modulus G given by the Hill averages 
(Table 1). It overlaps the original calibrated model result. 

martensite transformation during pseudoelastic loading to 
some macroscopic stress E. 

5.2.3. Discrepancy in open loop strain 

An important quantitative discrepancy is the underesti- 
mate of incremental strain per thermal cycle, AE cyc \ e , or 
“open loop strain” (Fig. 1). Although transformation is 
predicted to enhance plasticity, the simulations predict 
noticeable A E cycle values only at large X bias (>350 MPa) 
while experiments show noticeable values even at small 
Ebias (50 MPa). This discrepancy may be due to the nature 
of the “aggregate” constitutive relation (Section 3.2.3) and 
the assumption that the rate of slip activity (y,, Sec- 
tion 3.2.4) and rate of martensite formation (v,, Sec- 
tion 3.2.5) are both computed from the average aggregate 


stress. In reality, the local stress held around martensite 
plates can be sufficiently large to drive local plasticity, as 
evidenced by recent transmission electron microscopy of 
dislocation content in micron-scale, pseudoelastically com- 
pressed single crystals [27] and load-biased thermal cycling 
of larger single crystals [30], 

5.2.4. Assumptions concerning martensite elastic moduli and 
thermal expansion coefficients 

The simulations reveal several aspects concerning elastic 
moduli and thermal expansion coefficients. First, Fig. 8 
shows that the predicted transformation strain E T during 
thermal cycling is about 0.01 greater if a common assump- 
tion, C M ~ Vi C A , is used (e.g. [12,13]). This assumption is 
erroneous since both experiments [55,56] and first principle 
calculations [34,51] show that martensite is stiffer than aus- 
tenite (Table 1). A more compliant martensite increases E T 
by increasing the elastic strain and texturing in the mar- 
tensite; the former is deduced from the unloading path 
from A to A in Fig. 8. A reasonable approximation to 
the calibrated results is to use C M = C M(isotropic ), where 
C M (isotro P ic) is an isotropic matrix with Young’s and shear 
moduli given by Hill’s averages (see Table 1). The results 
are based on the anisotropic moduli of Hatcher et al. 
[34]. If values from Wagner and Windl [51] are used, smal- 
ler Hill's averages result (E h = 122 GPa, G h = 45 GPa), but 
Ca is also smaller. The net effect is a ~0.3% increase in E\ 
(Fig. 8), i.e. it is virtually unchanged. The observation that 
martensite can be represented by average isotropic elastic 
moduli suggests that it is also reasonable to use isotropic 
thermal expansion coefficients, at least for aggregate-scale 
simulations. This supports the adoption of isotropic 
descriptions for A, h . A and A th _ M in Section 4. 1 . 

5.2.5. Texture evolution 

The calibrated simulations are qualitatively consistent 
with measurements of textured martensite formation 



Bias Stress 2 bias (MPa) ^ Bias Stress 2 bias (MPa) 


Fig. 9. Normalized neutron diffraction intensities / vs. bias stress £ b i as for (a) (1 0 0) and (b) (Oil) martensite planes oriented parallel to the loading axis, 
for the 49.9 at.% Ni-Ti alloy. The intensities are measured at 0^ n = 30 °C, following thermal cycling to 0 ma![ = 230 °C, and are normalized by the intensity 
at Z bias = 100 MPa. The experimental data is from Ref. [2], The calibrated model parameters are presented in Table 1. Also shown is the calibrated model 
result using the interaction matrix h tu of Patoor et al. [23]. 
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during thermal cycling at X bias = — 150 MPa. During heat- 
ing to 1 30 °C, in situ neutron diffraction reveals peaks from 
retained martensite with plane normals N = [0 1 1] _L to the 
loading axis and also N = [1 0 0] || to the loading axis. Mar- 
tensite with N = [0 1 1] || or N = [1 0 0] _L is not observed 
[56]. The calibrated model predicts the same results at 
1 30 °C; in particular, there is 5% retained martensite, some 
of which has N = [0 1 1] _L and N = [1 0 0] ||, but none has 
N = [0 1 1] || or N = [1 0 0] X. This complements earlier 
work by Gao et al. [57], where isothermal loading was con- 
sidered instead of thermal cycling. 

Fig. 9 shows that the calibrated model qualitatively cap- 
tures the relative trends in [1 0 0] || and [0 1 1] || intensities 
with increasing X bias . The intensities are measured at 
0 m [ n = 30 °C and are normalized by values at 
2-bias — 100 MPa. The simulation results depend slightly 
on £p(pre) because the residual stress state aids or opposes 
the bias stress, depending on the location. For example, 
the normalized intensity [1 0 0] shifts downward by ~0.03 
for £p(pre) = 0 (vs. -0.7%) and [0 1 1] is unchanged. Results 
using the hardening matrix of Patoor et al. [23] (see the dis- 
cussion following Eq. (11)) tend to be less dependent on 
Ebias because the martensite texture is over-predicted at 
small X bias , leaving little room for additional texturing at 
larger £ bias . 

6. Conclusions 

Several key elements are required to capture the thermal 
cycling and isothermal deformation trends for polycrystal- 
line NiTi (49.9 at .% Ni). A critical feature is the coupling 
of transformation and plasticity. With appropriate data, 
the calibrated model captures (i) the effect of temperature 
on isothermal deformation loading response; (ii) the effect 
of bias stress on critical temperatures, transformation strain 
and martensite texture evolution during thermal cycling; and 
(iii) the effect of deformation processing in the austenitic 
state on the two-way effect and uniaxial loading at 0 > 0 AF . 

Two other key elements are the incorporation of aniso- 
tropic elasticity and an augmented form of the martensite 
hardening matrix h n , (Eq. (11)). The first corrects a fre- 
quent assumption that martensite moduli are those 
of austenite. The augmented h tu (Eq. (11)) increases the 
self-hardening term, thereby capturing the gradual increase 
in transformation strain and martensite texture with bias 
stress during thermal cycling. 

The simulations capture plasticity-transformation cou- 
pling in an aggregate sense but not on the martensite plate 
scale. On this smaller scale, retained martensite can affect 
texture evolution, and the stress redistribution in the 
vicinity of the plates may promote dislocation substructure 
[27-30] and affect the critical driving force for martensite 
formation (f c , Eq. (9)). Some consequences are discrepan- 
cies in the amount of reverse transformation during 
unloading and subsequent heat treatment (Fig. 3), the mag- 
nitude of ratcheting during thermal cycling under a bias 
stress (Fig. 1), and the effect of upper cycle temperature 


on transformation strain (Fig. 6) [2]. In principle, the use 
of 24 type II rather than all 192 predicted type I and type 
II habit plane variants [48] may affect martensite texture 
predictions, but it is unlikely to address the aforementioned 
discrepancies. 

Concerning future activities, the present work could ben- 
efit, in principle, from detailed studies of variant-variant and 
variant-plasticity interaction within aggregates. For exam- 
ple, the recursive method of Sittner and Novak [21,22] pro- 
vides a modeling-experiment approach with which to 
supplant the martensite hardening matrix (/;,„) in the present 
FE simulations. However, the present work also requires 
information about phase transformation-plasticity coupling 
at the subaggregate (discrete variant) scale. Currently, quan- 
titative information about such coupling is lacking. Inelastic 
deformation processes in martensite, such as detwinning and 
reorientation, can be addressed as described by Pan et al. 
[26]. In principle, martensite plasticity can also be incorpo- 
rated, but detailed experimental evidence of operative slip 
systems and critical stresses for activation is lacking. Given 
the recent observations of complex dislocation structures 
induced by transformation [25,29,30], it is remarkable that 
an aggregate-based, grain-scale model is so successful. 
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